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Abstract 


Principal curvatures and principal directions are exploited ubiquitously in shape analysis. On one hand, umbilical points severely 
hinder the analysis as the singularities in the vector field of principal directions; on the other hand, providing shape-intrinsic 
qualitative information about a surface, they are desirable quantities in some applications. Umbilical points are fundamental for 
geometric analysis, but their accurate computation on a discrete surface is still challenging. In this paper, we develop a simple and 
effective method to detect umbilical points on a triangle mesh, in which any parametrization works. In particular, we propose two 
practical processes for local parametrization by orthogonally projecting or conformally transforming the matrix onto a specified 
parametric plane. Furthermore, we make a systematic analysis on our method and demonstrate its convergence behavior. The 
algorithm of our approach is flexible and easy to implement for a triangular mesh of arbitrary topology. 
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1. Introduction 


Principal curvatures and their corresponding directions are 
basic entities for shape analysis because they entirely exhibit 
the variation of curvature across a smooth two-manifold sur- 
face in a three-dimensional Euclidean space. Their calculation 
is a fundamental issue for enormous geometric applications in 
many aspects such as shape retrieval, shape registration, surface 
remeshing, surface segmentation, and shape design. 

An umbilical point is a singularity with anomalous behav- 
ior in the vector field of principal directions. It is of theoretical 
interest because lines of curvature, which integrate the field, be- 
come complicated near it. In applications favoring the orthogo- 
nality of the field, such as feature-aligned remeshing [1, 9], um- 
bilical points make waves; we have to isolate them by design so 
as to deal with them separately. Umbilical points are invariant 
features of the surface’s differential structure, acting like finger- 
prints across the surface [10]; hence, their reliable determina- 
tion is also significant for quantitative analysis in other settings, 
such as for describing or identifying objects [16], shape inter- 
rogation [13], ship-hull design and progressive additional lens 
design [17]. Although a wealth of estimators for differential 
quantities have been proposed in the vast literature of discrete 
geometry, the literature lacks a critical examination regarding 
how to detect umbilical points over a mesh as opposed to a s- 
mooth surface. 

The purpose of this work is to investigate a mathematically 
sound and computationally efficient approach to umbilical point 
detection over a mesh model of arbitrary topology. The remain- 
der of this paper is organized as follows. In the next section, a 
summary of related work and their mathematical background 
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are provided on the numerical computation of umbilical points 
over a mesh. In Section 3 mathematical preliminaries in classi- 
cal differential geometry are described concerning the charac- 
teristics of umbilical points and shape operators; and then the 
theoretical foundation of our method is developed to linearly 
interpolate the shape operators for umbilical points. We pro- 
pose our algorithms to re-formulate shape operators by a local 
transformation onto a parameter plane to facilitate the linear in- 
terpolation of shape operators in Section 4. Section 5 presents 
a unified error analysis on our method. The last two section- 
s present some experimental results with discussion and then 
conclude the paper, respectively. 


2. Previous work 


Umbilical points have been studied intensively. We shall re- 
state those work from the view point of numerical computation. 
Interested readers can refer to [5, 12, 13, 16] and references 
therein for more details. 

In the following text, a bold uppercase letter represents a 
matrix and a bold lowercase letter represents a column vector. 


2.1. An umbilical point 


Given a unit normal vector at a point p on a smooth sur- 
face, a normal plane at p is one that contains the normal and 
therefore also contains a unique direction tangent to the sur- 
face, cutting the surface in a plane curve. This curve usually 
has different curvature for different normal planes at p. The t- 
wo principal curvatures at p, denoted as x; and x2, are the max- 
imal and minimal values, respectively. The pioneering Swiss 
mathematician, Leonhard Euler, first discovered that the two 
principal directions of a; and q@2 at p, corresponding to the tan- 
gent directions realizing x; and x2, respectively, are mutually 
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perpendicular. The normal curvature «x, in the tangent direction 
making an angle 0 with a, is characterized by the well-known 
Euler’s formula: 


Ky, (0) = ky cos? 0 + K2 sin? @ 


If kı = ko, however, K, in all directions is the same; this point is 
called an umbilical point and each direction is principal. Um- 
bilical points exhibit special properties and their computation 
is an important topic in shape analysis. Two commonly used 
approaches are related to two properties of umbilical points, re- 
spectively: 


1. In classical differential geometry, an umbilical point is de- 
fined as a point where x; and x are equal. A trivial algo- 
rithm to locate an umbilical point is to locally minimize 
f = |Kı — k2| at a vertex against other vertices adjacent to 
it, or, to have the ratio of x; and x2 larger than a thresh- 
old [9]. Its equivalent form of g = i — Kg is more fre- 
quently exploited in practice, where km = (Kı + K2)/2 is 
the mean curvature and Kg = k,K2 the Gaussian curvature. 
Because the notion of continuous geometric attributes on a 
smooth surface do not match those on its piecewise linear 
counterpart (more or less), it is uncommon to determine an 
umbilical point by f or g alone on a sampled surface. This 
method is more widely utilized for an analytical surface on 
which numerical optimization can be performed to situate 
an umbilical point accurately. For instance, a method was 
proposed in [12, 13] for a free-form surface by solving a 
system of polynomial equations. For an implicit surface, 
novel evaluation functions and numerical algorithms have 
been developed to extract umbilical points in [20]. 

2. Another common method is based on the index of prin- 
cipal direction field at a specified point p. The index of 
an umbilical point describes the way the lines of curvature 
turn around the umbilical point and can be computed by: 
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where @ is the angle through which the coordinate axes 
would have to be rotated to align them with the local prin- 
cipal axes at p [2]. The integral is taken over a small coun- 
terclockwise circuit c around p. In implementation, the 
computation is simplified by pulling the field back onto the 
tangent plane at p and computing the circuit there [13, 16]: 


1 
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where @ is the angle between a direction of the field and 
some fixed direction in the parametric plane. The integral 
is now taken going around the circuit c around p but in the 
plane. 
The index of a generic umbilical point is either of index 
-3 or of index + L, Thus this method provides a conclusive 
test for an umbilical point. If the index is zero no umbilical 


point is present. 


2.2. Linear interpolation for an umbilical point 


Localizing an umbilical point is more than just a matter of 
evaluating |k; — k2| or calculating the index. Because it is likely 
that the actual umbilical point will fall between the sampling 
points, the structure of the principal direction field should be 
taken into account [16]. A method by linearly interpolating 
some quantities at vertices of a geometric primitive — for ex- 
ample, a triangle over a mesh — is being sought to yield an 
umbilical point. 

In a2D symmetric tensor field T of type 2x2, the location of 
a singular point on a discrete grid is reduced to an interpolating 
point of zero using bilinear interpolation of the tensor-matrix 
components between the grid vertices [6]; later generalized to 
linear interpolation within each triangle Aabe by solving a lin- 
ear system of T [18]: 


xD, + yDp + (1 -x-y)D,. = 0 (3) 


where the deviator part D of T is defined as D = T — Str(T)E; 
here E is an identity matrix of 2 x 2. 

This method is then employed in remeshing for an umbilical 
point over a triangular mesh by flattening it to a 2D parametric 
plane and by interpolating the curvature tensor [1]. For a para- 
metric surface, linear interpolation of curvature tensors in the 
parametric domain is advantageous to fixing an umbilical point 
inside a primitive cell rather than only on a vertex as the former 
two do. Their practice, however, has to rely on a globally con- 
formal parameterization of the mesh in order for symmetrical 
curvature tensors in the parametric domain. Parameterization is 
a central issue in computer graphics and it is a nontrivial task to 
obtain a conformal one. Furthermore, a parametrization is not 
indispensable in most applications. Needless to say, the global 
parametrization requires the mesh surface to be equal to a pla- 
nar disk or a planar square in topology; additional errors will be 
produced in the parametrization, too. 

How to estimate an umbilical point with a local parameter- 
ization is the key to popularizing the method through linearly 
interpolating certain quantities at vertices of a geometric prim- 
itive. The contributions in this work is such a novel method 
of local computation for umbilical points on a triangular mesh. 
It significantly frees the interpolation from both symmetrical 
tensors and conformal parameterization so that the underlying 
surface is of arbitrary topology. We also carry out a mathemat- 
ical analysis on it by examining its numerical and convergence 
behaviors of estimating an umbilical point on a discrete mesh 
surface — a method that has not yet appeared in the extant lit- 
erature. 


3. Mathematical preliminaries 


Although our method can be portrayed compactly, its math- 
ematics does not described explicitly in the literature. To pro- 
vide a practical foundation, a brief review of some important 
mathematical preliminaries related to the algorithm is in order. 


3.1. Background of differential geometry 

We summarize the relevant definitions and results about lines 
of curvature on a parametric surface to depict our method. Dif- 
ferential geometry in parametric representation has been devel- 
oped richly, and most facts can be found in a standard textbook 
on classical differential geometry when proofs are not given ex- 
plicitly. 

Given a variable vector u = [u,v] , in which u and v are lo- 
cal curvilinear coordinates, a vector-valued parameterized sur- 
face in terms of u is represented as 


Jf 


r(u) = [x(u, v), y(u, v), z(u,v) |" 


If r is regular, r, and r, are linearly independent and form a 
basis of the tangent plane at r(u). The shape operator S for r(u) 
is: 
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damental forms in the basis r, and r,, respectively; n is the unit 
normal vector. 

A basic property of the shape operator is Sa; = Kaj, i = 
1,2, where x; and œ; are the principal curvatures and the corre- 
sponding principal directions. 

The differential equation of lines of curvature in coordinates 
is given by 

Udu? + Vdudv + Wdv* = 0 (4) 
where U = EM-FL,V = EN-GL, and W = FN-GM. 


Therefore, we have a well-known result of an umbilical 
point in classical differential geometry: 


Lemma 3.1. r(u) is an umbilical point if and only if U = V = 
W =0Oatu. 


It seems that U = V = W = O are over-determined to situate 
an umbilical point; nevertheless, they are not independent, as 
shown by the following corollary: 


Corollary 3.2. r(u) is an umbilical point if and only 7 a 
at u. 

V=0 =X 
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Remark 3.3. In the sketch proof of Corollary 3.2, we use a 
concentrated notion of & = 4 = X for the sake of a compact 
proof. It infers a presupposition that M must be zero if F is zero 
for this equation to hold. A reader can examine that it does not 
affect the whole validity of the proof. We should also note that 
both E and G are always positive numbers for a regular r. 


Corollary 3.2 is inconvenient to characterize an umbilical 
point in numerical computation for a discrete mesh model. In- 
stead, defining t = [S11 — S22, S12 + Sal’, a more intuitive and 
practical condition can be derived for this purpose, say: 


Lemma 3.4. p is an umbilical point if and only if t = 0. 
Si; -S22 =0 V=0 We 


Si +S; =0 U-W=0 ` 
verify this lemma immediately following from Corollary 3.2. 


Proor. t = 0 e 


To the best of our knowledge, Lemma 3.4 is not described 
explicitly in the literature. Its power will be reflected by our 
new algorithm to compute umbilical points over a two-manifold 
mesh model with arbitrary topology. 


3.2. Linear interpolation of shape operators for umbilical points 


In this subsection, we further give an essential link that ex- 
ists between a shape operator and a curvature tensor under a 
conformal parameterization. We then develop a novel approach 
to linear interpolation based on shape operators for umbilical 
points rather than curvature tensors, which can relax the restric- 
tion of the symmetry of shape operators as opposed to curvature 
tensors. 


3.2.1. A shape operator and a curvature tensor 

In contrast to a shape operator, x; and œ; (i = 1,2) together 
are also thought of as the curvature tensor C = kiaia? +K20207 
in the computer graphics community. The tensor field T in 
Eq.(3) is set as C for umbilical points, but C must be symmetri- 
cal, which is too strict to exploit in applications. Because both 
S and C are closely related to x; and a;, i = 1,2, there must 
be a connection between them and a well-founded thought is to 
replace C with S so as to take full advantage of shape operators, 
which have been well studied in differential geometry. 

The following proposition reveals the relationship between 
S and C: 


Proposition 3.5. C = S if both are represented in conformal 
curvilinear coordinates. 


Proor. Let (u,v) be conformal coordinates. a; and a2 are unit 
vectors of principal directions in coordinates. Because a, -a2 = 
0, qia? + aa} = E holds. Thus C = kaja? + kana} = 
S (e107 + ana) =S. 


U - W = 0.2.2. A new method for umbilical points 


Lemma 3.4 and Proposition 3.5 imply not only that Eq.(3) 
can be simplified by t, but that it can also be applied to any 
parameterization so as to free S from symmetry, namely there 
being (x, y) € [0, 1] x [0, 1] and x+y < 1 fora triangle of Aabe 
such that 


xta + ytp + (1 — x- yte = 0 (5) 


We call Aabe meeting Eq.(5) as a zero-crossing triangle. 
A triangle containing an umbilical point is called an umbilical 


triangle. A zero-crossing triangle is unnecessarily an umbilical 
triangle in practical compuation. 

We shall continue the numerical analysis on Eq.(5) and il- 
lustrate more of its behaviors by experiments later — in Sec- 
tions 5 and 7, respectively. 


4. Finding a zero-crossing triangle 


Given a mesh model, a great many methods have been de- 
veloped to derive the principal curvatures and directions at each 
vertex [7, 19]. In this section, we develop a local parameteri- 
zation of high efficiency to synthesize the estimated principal 
curvatures and principal directions to be shape operators con- 
sistently in a coordinate system. 


4.1. Overview of the method 


Suppose principal curvatures and principal directions, i.e. 
shape operators at their respective frames, at vertices have been 
acquired by some means, either analytically or numerically. To 
depict our method, it suffices to carry out Eq.(5) on one triangle 
Aabc of the mesh model. Now imagine a right-handed uvw 
coordinate system at a point po inside Aabe with w-axis parallel 
to the normal to Aabe. We have a demand to project shape 
operators of a, b and c onto the uv plane, respectively. Once 
they are determined, interpolation is performed via Eq.(5). The 
resultant interpolating umbilical point is calculated as xa + yb + 
(1-x-vy)e. 

The rest of the issues involves projecting relevant shape op- 
erators to the parametric domain on Aabc. We shall portray our 
method to parameterize a shape operator locally for this pur- 
pose in the next subsection. 


4.2. Local parameterization of S 


We here illustrate the concrete implementations of our method 
with two local parameterizations. 


4.2.1. Local orthogonal projection 

We choose uvw-coordinates so that a point on the surface is 
at the origin and so that the three axes are the principal direc- 
tions and normal, respectively. The surface is then expressible 
in Monge form as 


1 1 
w= z% u? +k Vv") + zeo + 3biu’v+ 3bzuv? +b3v°) + O(u)*(6) 


in which O(u)* is a Peano remainder. 

For each vertex of the mesh, the Monge form in its neigh- 
borhood can be established numerically and uniquely to degree 
2 by estimating its normal, principal directions and principal 
curvatures from a set of near-neighbor vertices. However, the 
shape operators of two distinct vertices cannot be put into arith- 
metical operation directly because they are obtained separately 
and measured in their respective frames. To reconcile them to 
a consistent frame, the Monge forms need to be reshaped in a 
common parametric domain. Given a local xyz coordinate sys- 
tem at po and a Monge form in the local uvw coordinate system 
at px as shown in Fig. 1, we intend to re-parameterize S|p, from 


Figure 1: Coordinate transformation for local parameterization of S. Frame p- 
uvw is a Monge frame at px to be re-parameterized. P% is the projection of px 
on the xy plane. 


in the frame p,-uvw to in the frame po-xyz. As described in Sec- 
tion 3.1, S is computable by setting r(x) = [x, y, z]” if Jx(z) and 
H(z) are computable, where Jx(z) and Hx(z) are the Jacobian 
and the Hessian matrices of z in terms of x, respectively. 

Let ei, e2, and e3 be the unit vectors in the frame po-xyz 
along the positive directions of the x, y, and z axes, respectively; 
€,, &5, and ë; are those of the u, v, and w axes, respectively, as 
displayed in Fig. 1. Let aj; = e;-€;, bi = (pe—po)-ei, i, j = 1,2, 3, 
and U = [u,v,w]’, X = [x,y, zl’. We have 


X = [aij], U + Bila (7) 


Let a = [ay2,a11]", b = [a2,a2]", € = [a23,a13]’, A = 


-1 


a a Kk, O 
[b, —a], B = | py and K = | a © | By elemen- 
tary computation, we have 
Ju) = B (8) 
H,(u) = det([c,b])det(B) ATKA (9) 
H,(v) = det({a,c]) det(BXATKA (10) 
Jw) = 0 (11) 
H,(w) = BKB (12) 


where = means the evaluation is taken at p (or p,), e.g. 
J,Wlp: =0. 

According to z = a3;U + a32v + a33w in Eq.(7) and the resul- 
tant Eqs.(8)-(12), we eventually achieve J,(z) and Ay(z) at px; 
and thus Sj, . 


4.2.2. Local conformal transformation 

Instead of orthogonally projecting them onto a parametric 
plane, we can flatten the matrices of shape operators onto the 
triangle plane by transforming the vertex frames to the face 
frame rotating around cross product of their respective normal- 
s with the plane normal, e.g e3 and €; in Fig.l, while keeping 
the angles of the principal directions [15]. Only three trian- 
gle vertices account for the interpolation and this local confor- 
mal transformation can regarded as a local structure of a certain 
global conformal mapping. Our experiments show that it works 
almost as good as the local orthogonal projection. 


5. Error Analysis 


We now set out to address the fundamental question of its 
convergence rate as the mesh resolution is increasing. 


5.1. Invariant in linear interpolation 


We shall examine the quantitative behavior of our linear in- 
terpolation method on a surface locally represented in the Mon- 
ge form. Before proceeding, we give a lemma to this end: 


Lemma 5.1. The interpolating point of a zero-crossing trian- 
gle is invariant in any curvilinear coordinate arising out of a 
rigid transformation. 


Proor. S is expressed in terms of u. Given a rigid transforma- 
cos@ —sin@ 
sin@  cosé 
and Ug is a constant translation, then S in terms of û can be rep- 
resented by R’SR. The elementary computation shows that 


tion between u and û, i.e. u = Rû+uọ, where R = | 


tû) = Mtu) (13) 
cos? 0 — sin?  2cos@sind 
—2cos@sin@ cos? 0 — sin? 0 
Suppose the three vertices of a zero-crossing triangle are 
represented as (a, b, c) in terms of u and (â, D, ĉ) in terms of û, 
respectively. We have 


where M = 


xt(a) + yt(b) + (1 — x- y)t(®@) = M (xt(a) + yt(b) + (1 — x - y)t(o) 


Because det (M) = 1, it holds that xt(a) + yt(b) + (1 — x - 
y)t(€) = 0 © xt(a) + yt(b) + (1 — x — y)t(c) = 0. We verify this 
lemma immediately. 


5.2. Error Analysis 


By substituting the Monge form of Eq.(6), elementary com- 
putation works out as follows: 


Si Si 
S = 
| Sa S2 | 
_ K1 + bou + biv biu + bov 2 
= biu + bov Ko + bou + b3v + Olu) 
(14) 
Hence 
_ | Ki -K+ (bo — bo)u + (bı — b3)v 
M S | (biu + bv) | 
= ty +ut; +vt, + 0? (15) 
where tọ = | a2 |t = | a b2 and tz = | a bs | 
1 2 


tı and tz are the principal linear parts. 

Given a Aabc in the uv domain with origin p, it covers 
a point of d = [ū,ọ]7 under consideration. Without loss of 
generality, we can imagine that p is inside Aabc, for example, 


the centroid, according to Lemma 5.1. Suppose the barycen- 
tric coordinates of d in terms of Aabe were (x,y). Denote 
h = max(\a — pl, |b — pl, |c — pl), then 


xt(a) + yt(b) + (1 — x- y)t(c) = to + at; + vtz + O(h?)(16) 
According to Eqs.(15) and (16), we have 


t(d) = xt(a) + yt(b) + (1 — x—y)t(c) + O(n’) (17) 


Eq.(17) reveals the operational principle of our method. It 
means that t is approximated by linear interpolation to accuracy 
of O(h). Here h can be relaxed to be the resolution size of the 
mesh — the maximal length of all edges. It also reflects that 
the key to linear interpolation is that t can be linearized locally 
such that t is dominated by its linear terms. Either tı + 0 or 
t2 + 0 is required to make Eq.(17) in proper operation, namely 


bıb2 #0 or (bo — b2)(b1 — b3) #0 (18) 


Eq.(18) is satisfied with a generic umbilical point. If the 
mesh is fine enough, a zero-crossing triangle of Aabe implies 
that Aabc is also umbilical by both Eqs.(17) and (18). 

Because an interpolating point can approximate a zero point 
of t — that is, an umbilical point — we need to address the 
further issue of what the accuracy is between an umbilical point 
and an interpolating point to complete our discussion. 


Assume d be an interpolating point, we have by Eq.(16): 
to + it) +t = O?) (19) 


Suppose t(d) = 0 inside Aabe — that is, d = [u,v]? is an 
umbilical point — then 


to + ut; + vty = O(n’) 
Eq.(20) subtracted from Eq.(19) gives 
(a — u)ti + (6 — v)t2 = O(N’) 


(20) 


According to the Cauchy-Schwarz inequality, we have 


la -ā| = Vā - u)? + © - v}? > Oh’) 


On the other hand, if the umbilical point is generic and the 
mesh resolution is high enough, t is approximated to sufficient 
accuracy by its principal linear components and thus a zero- 
crossing triangle coincides with an umbilical triangle. In this 
case we obviously have ja — d| < O(h). Therefore, we have the 
following proposition. 


Proposition 5.2. If the mesh resolution is sufficiently high, the 
interpolating point of a zero-crossing triangle approximates an 
umbilical point to accuracy between O(h) and O(h?). 


5.3. Error propagation 

The above analysis is based on the exact S. However, prin- 
cipal curvatures and principal directions are estimated on accu- 
racy in practice up to some extent. Thus, their errors have to 
be considered in estimation of umbilical points by linearly in- 
terpolating S, i.e. error propagation from shape operators to an 
interpolating point. According to the discussion in Section 5.2, 
it entirely relies on the accuracy of the estimated S. 

We think of it under the conditions presented by [8]: 


Proposition 5.3. Given the position, gradient, and Hessian of 
coordinate functions that approximated to O(h“*!), O(h®), and 
O(h*-'), respectively, and assuming the condition number of the 
Jacobian matrix is bounded, we have: (a) The angle between 
the computed and exact normals is O(h*); (b) the components 
of the shape operator and curvature tensor are approximated to 


o(ht"!). 


For a mesh model, d is typically 2 such that S is approxi- 
mated to accuracy O(h), which gives us the following corollary: 


Corollary 5.4. A resultant umbilical point by linear interpola- 
tion is typically approximated to accuracy O(h). 


6. Recondition umbilical points 


Umbilical points are difficult to be reconditioned with com- 
plete fidelity although they are inherent in a surface. There 
are two factors of S tensor distribution as well as geometry de- 
tail across the strongly affecting their computational quality. In 
practical usage, two crucial issues are to be addressed: One is 
how to judge an umbilical point to be mistakable and the other 
is how to compute operable umbilical point set. 

Judging mistakable umbilical points. To the best of our 
knowledge, no work has been done for this purpose in the lit- 
erature yet. We propose a method to determine if an umbili- 
cal point is mistakable by a necessary condition: An umbilical 
point is either spherical or planar and thus distributes over an 
elliptic (kg > 0) or a planar (kg = 0 and ky = 0) regions of the 
surface [5]. It generally occur as an isolated point in the ellip- 
tical region and, furthermore, a flat region is left out of account 
since it is unusually of interest in actual applications. Another 
obvious observation indicates that a significant umbilical point 

should be supported over a sufficiently large elliptic region; in 
- contrast, the local shape variation only yields a small region of 
same type. Based on those two considerations, our first strategy 
is devised as follows: 


1. First, the surface is partitioned into parts with surface type 
labels based on xg and xy. Only elliptic surface types are 
remained, i.e. Peak and Pit types! [4], for further treatmen- 
t. 

2. Then, two fundamental operations in mathematical mor- 
phology — dilation operator D and erosion operator E, are 
used in turn. They compose an opening operator D” o E” 
with n erosions followed by n dilations and applied to the 
elliptic parts to erase narrow connectors and open up large 
holes or caves. n is determined by the the sizes of narrow 
connectors to be erased. Umbilical points are eroded away 
near the boundaries of elliptic parts. 

3. Finally, the remaining umbilical points are kept over the 
elliptic parts. 


Inconsequential umbilical points can be clean out significantly 
in this way. We will demonstrate this effect in our experiments. 

Computing operable umbilical points. A real-world mesh 
model usually contains a plenty of curvature variation across its 
surface arising from various desired or undesired factors, such 


as detail enhancement and sampling noise. The sensitivity of 
umbilical points to those variation reflects their characteristic 
response to the shape undulation, leading to a massive amount 
of umbilical points. But too many umbilical points might af- 
fect their application. To make them operable in practice, their 
number should be drastically reduced. For example, it is re- 
duced to simplify the topology of the extracted field so as to cut 
a surface for quadrangulation [3]. In this scenario, an approach 
to smoothing or denoising shape operators will be beneficial to 
minimizing the number to meet different requirements — such 
as removing superfluous points by C smoothing [1] or by refine- 
ment of principal frame fields [16]. In contrast to those methods 
by denoising 2-order differential quantities, we smooth normals 
by a linear filter and shape operators are then estimated based 
the resultant normals. Our experimental results will show that 
our method can decimate undesire umbilical points effectively. 


7. Experimental results and discussion 


In this section, we shall demonstrate our method by experi- 
ments and comparisons with other methods based on a surface 
of analytical representation and some real-world mesh models. 


7.1. Method validation 


We confirm our proposed techniques by using a practical 
example with known umbilical points on a given Monge patch: 


z = 0.5(x? +y) + 1.257 - 0.5x°y + 0.2xy? + 1.6? (21) 


Since Eq.(21) is a polynomial expression, all its umbilical points 
can be found by a solver of polynomial system, as done in [20]. 
As a matter of fact, there are exactly three umbilical points on 
the patch: 


pı = [0,0,0]" 
p2 = [-0.95695088, —0.18571175, —0.50829269]" 
p3 = [-0.39165274, —0.81922023, —0.52924773]" 


Given a set of 5.4K points by a random sampling inside 
the region [—1.2, 1.2] x [—1.2, 1.2] on the parametric plane, a 
Delaunay triangulation is constructed. We polygonize such a 
patch to be a triangular mesh using the Delaunay triangulation. 
For the sake of clarity, pı, p2 and p3 are marked noticeably by 
their respective bounding boxes as in Fig. 2. 

We devise the experiments following the line from fine shape 
operators to coarse ones. 


7.1.1, Numerical behavior 
We first consider an instance excluding the influence of noise. 

The exact principal curvatures and principal directions are both 
computed analytically at each vertex according to Eq.(21) so 
that shape operators are ready to be applied to our approach re- 
gardless of numerical errors in discrete estimation. The linearly 
interpolating points are displayed in Fig.4. Not unexpectedly, 
it exhibits favorable interpolation despite the irregular mesh: 
three umbilical triangles are figured out correctly and the inter- 
polating points are close to the ground truth umbilical points for 
sufficient accuracy. 


Figure 2: A designed irregular mesh of a surface patch on the parametric 
domain with ground-truth umbilical points of pı, p2 and p3, marked inside 
their respective boxes with same indexes. The zoom-in views are exhibited in 
Fig.4(a)-4(c) of umbilical points resulting form index of vector field, curvature- 
tensor interpolation and shape-operator interpolation, where shape operators 
are computed analytically. Those with shape operators estimated in discrete 
fashion are exhibited in Fig.5(a)-5(d), respectively; umbilical points inside Box 
4 are suspicious. 


Figure 3: Sharply changing principal directions near erroneous umbilical points 
in the parameter plane by local orthogonal projection. 


We then compute the shape operators by a discrete estima- 
tor: We approximate the rough normal at a vertex using the 
method of [14] and choose it as the z axis to fit a polynomial 
surface in its two-ring neighbors to estimate the applied nor- 
mal; based on these resultant normals, the shape operators are 
finally computed using normal-fitting method. Fig. 5 exhibit- 
s the results using the discrete estimator of shape operator. In 
contrast to exact S, the resultant p; lies slightly outside the rel- 
evant umbilical triangle and thus the zero-crossing triangle of 
pı fails to coincide with its umbilical triangle (See Fig. 5(a) for 
zoom-in views). But it is acceptable since the ground truth pı 
is near a triangle edge. Moreover, there are some exceptions 
for interpolation of discretely estimated shape operators: Two 
additional umbilical points are produced in error, as shown in 
Fig. 5(d). They obviously come into conflict with the ground 
truth. We elucidate the reason for these wrong points. We 
project the principal directions onto the parametric plane, as 
illustrated with Fig. 3, and we can see that principal directions 
change very sharply near them. It is really an ambiguous case 
due to the sudden change of the principal directions. The el- 
ementary calculation shows that the angle between the nearby 
principal directions near it is right-angled as 92°. They lie at the 
transition region of different surface types and can be filtered by 


an opening operator. 


Table 1: Linear interpolation with ground truth shape operators based local 
orthogonal projection(O) and local conformal transformation(C). 


Estimated Projection 
1.10570e-3,7.25923e-4,8.05220e-4 O 
Pl | 710569e-3,7.25909e-4,8.05220e-4 C 
-0.959762,-0.184829,-0.515272 O 
P2 | 0.959762,-0.184829,-0.515272 C 
-0.388567,-0.827572,-0.551645 O 
Ps | -0.388567,-0.827572,-0.551644 C 


Table 2: Linear interpolation with estimated shape operators based local or- 
thogonal projection(O) and local conformal transformation(C). 


Estimated Projection 
2.90345e-3,2.18191e-3,8.30664e-4 (0) 

Pi [2.90349e-3,2.18194e-3,8.30665e-4 C 
-0.962192,-0.185054,-0.520515 (0) 

P2 | 0.962192,-0.185054,-0.520515 C 
-0.388707,-0.827926,-0.552492 (0) 

p3 | -0.388707,-0.827926,-0.552492 C 


We have shown that any valid parametrization can apply to 
the shape-operator-based interpolation. Besides local orthog- 
onal projection, we also give a local conformal transformation 
as a means of implementing the method in Section 4.2.2. Ta- 
ble 1 and 2 list the experimental data of the two implementa- 
tions. There is no appreciable difference between them, also as 
expected in our more experiments. Therefore, shape-operator- 
based interpolation via orthogonal projection is only taken in 
following comparison experiments. 


7.1.2. Comparisons 

Can other comparable methods do better than our proposed 
method since ours might fail to respect the true umbilical points 
as shown in Fig. 5(d)? In this subsection, we compare our 
shape-operator-based interpolation method to the index-based 
method grounding on Eq.(2) and the curvature-tensor-based in- 
terpolation method via global conformal parameterization, de- 
noted with S, I and C methods for short, respectively. For 
curvature-tensor interpolation, we apply LSCM parameteriza- 
tion in [11] to performing the global parameterization; because 
it requires the mesh to be homeomorphic to a planar disk or 
a planar square in the parametric domain, the Monge patch of 
Eq.(21) happens to meet. In this scenario, we can re-flatten it to 
the planar domain of conformal parameterization. 

Reusing the shape operators in Section 7.1.1, both analyti- 
cal and discrete, we compute the umbilical points with the three 
methods, respectively. Figs.2 and 4 illustrate the resultant inter- 
polating points in analytical case. So do Figs.2 and 5 in dis- 
crete case. It is obvious that the interpolation behavior of C 
and S methods are quite similar — something that is expect- 
ed. Without regard to the error effects of global discrete pa- 
rameterization, their resultant umbilical points can be thought 
of comparable; unfortunately, the error of approximately global 
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Figure 4: Resultant umbilical points based on exact shape operators near p1, p2 and p3. A green °’ represents an exact umbilical point, a magenta °x’ results from 
linearly interpolating curvature tensor, and a blue ’+’ is from linearly interpolating shape operators; a red ’o’ is determined by index-based method. 
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Figure 5: Resultant umbilical points via index of vector field, curvature-tensor interpolation with global conformal parameterization, shape-operator interpolation 
with local parameterization, based on discrete estimate. 


parametrization can make the interpolation less reliable, which 
will be shown by experiments in the following text. In vivid 
contrast against both C and S, I method can yields a cluster of 
points near an umbilical points rather than one. 

Fig.5 makes a response to the inquiry at the beginning of 
this subsection: both I and C methods cannot do better than S. 
These facts will be more evident in the following experiments. 


7.2. Linear interpolation over mesh models 


We shall demonstrate our method with real-world mesh mod- 
els. 

Fig. 6 shows the resultant umbilical points on a hand mod- 
el by the three methods, respectively. To help understanding 
their distribution over the mesh, one branch of estimated prin- 
cipal directions are also displayed. In the I method, a cluster of 
vertices can be detected arising out of an umbilical point, typ- 
ically either 1, or 2, or 3 (See Fig. 6(a)); the C method fails to 
work near the tips of the fingers except for the thumb because 
the parameter gradients of the underlying global conformal pa- 
rameterization are too large (See Fig. 6(b)); the S behaves the 
best(See Fig. 6(c)). And, of course, S method will work better 
if the model is decomposed into more smaller charts and then 
flatten to reduce the distortion. 

The second example is of a rocker-arm model. It is close 
and thus the global parameterization is not performed. Fig. 7 il- 
lustrates the difference between the S method and the I method. 
The former usually determines only one umbilical point near 
a singular point in the principal direction field, but the I does 


(b) A zoom-in view 


Figure 7: The rocker-arm model. The detected umbilical points by the I method 
are drawn as a smaller red ball while those by the S are larger black balls. 
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Figure 6: The hand model. Local methods, i.e. I and S, behave themselves, but the global C results in some fallacious umbilical points near the tips of the four 
fingers other than the thumb. (d) shows the resultant points processed by elliptic regions based on S. 


not. Because sufficient points per boundary curve are required 
to estimate the index, the index-based is more sensitive to the 
irregularity of the mesh. 


7.3. Filtering umbilical points 


Given trivially-estimated shape operators, the previous ex- 
ample with ground-truth umbilical points exhibit some mistak- 
able ones even if the mesh is noise-free(See Figs. 2 and 5(d)). 
As pointed out in Section 6, two strategies can be exploited to 
recondition the resultant umbilical points and thus to increase 
the practical utility of those resultant points. 

The first one is applying the opening operator to an ellipti- 
cal region over the surface so as to erase narrow regions usual- 
ly related to geometry details or transition regions of different 
surface types. The bigger n is, the more umbilical points are 
removed. To balance the fidelity and practicability of resultant 
umbilical points, n = 2 is chosen in our experiments. Further- 
more, a flat region is left out of account since it is unusually 
of interest in practical applications. We apply this strategy to 
Fig. 5(d), and the result is shown in Fig. 8. We can see that 
the mistakable points outside elliptical regions are eliminated. 
Fig. 6(d) shows the reconditioned umbilical points for the hand 
model. 

However, the first strategy is not capable of handling com- 
plex geometry, for instance, the Standford Bunny with plenty 
of curvature variation. Fig. 9 shows the interpolating umbili- 
cal points handled without and with opening operator, respec- 
tively. Because the ’low’ mesh resolution of the Bunny cannot 
accommodate its surface roughness, too many umbilical points 
still remain even after opening operation, which will affect their 
practical utility. In this setting, an approach to smoothing or de- 
noising shape operators will be beneficial to make further im- 
provement on the filtering quality. As an alternative, we apply 
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Figure 8: Recondition the umbilical points on the surface patch by surface type. 
(Refer to Figs. 2 and 5). Colorful regions are elliptical regions where umbil- 
ical points should locate. Mistakable umbilical points in Box 4 are excluded 
successfully. 


linear fitting to the normals before shape operators are estimat- 
ed. The experimental results are shown in Fig. 10. This strate- 
gy works very well for the rocker-arm model; as illustrated with 
Fig. 11, it makes a significant improvement in contrast to Fig. 7. 

Fig. 12 results from the combination of the two strategies. 
It shows that the resultant points reflect symmetry for a sym- 
metrical shape. They can be used as anchor points to extract 
the symmetrical plane of the Torso. 


7.4. Computational cost 


Given principal curvatures and principal directions of ver- 
tices over a mesh, linear interpolation of shape operators for 
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Figure 9: Umbilical points over the Bunny. The top calculates umbilical points 
straightforwardly; the bottom filters the points by elliptic regions. 
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Figure 10: Umbilical points over the Bunny by smoothing shape operators with 
normal fitting. 


Figure 11: Umbilical points over the rocker-arm. Umbilical points over the flat 
regions are discarded. 


Figure 12: Umbilical points over the symmetrical Human-Torso model. The 
resultant points reflect the symmetry of the model. 


umbilical points is performed by going over all triangles. For 
an interpolation of each triangle, no optimization computation 
is required. Instead, only some necessary floating point arith- 
metic for local reparameterization of shape operators onto each 
triangle plane (as described in Section 4.2) is needed. There- 
fore, the computational cost of our algorithm increases linearly 
with the number of triangles and runs on-the-fly as opposed to 
the cost for estimating shape operators. 


8. Conclusion 


In this paper, we described a novel approach developed for 
umbilical-point detection over a mesh by interpolating shape 
operators via a local parameterization. The shape operators do 
not need to be symmetrical any more and we propose using a 
local parameterization by orthogonal projection or local confor- 
mal transformation to cooperate with the interpolation of shape 
operators. Furthermore, we conducted a systematic error anal- 
ysis on the interpolation to demonstrate its efficiency. 

Linear interpolation is superior to the index-based method 
because the latter should collect sufficient neighboring points 
for index computation and, accordingly, the former is less sen- 
sitive to the regularity of the underlying mesh. In contrast to 
curvature-tensor based on global conformal parametrization, the 
proposed method in this work is significantly advanced because 
its local parametrization is done more robustly and no addi- 
tional errors are introduced in parameterization; furthermore, 
it does not suffer from the topological requirements of the mesh 
surface for global parametrization; and thus significant com- 
putational burden and implementation complexity are removed 
without the construction of a global parametrization or the seg- 
mentation of a complex model. As a result, our method saves 
considerable time, and is more flexible and easy to implement. 
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